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We develop a steady-state analytical and numerical model of the optical response of power-recycled 
Fabry-Perot Michelson laser gravitational- wave detectors to thermal focusing in optical substrates. 
We assume that the thermal distortions are small enough that we can represent the unperturbed 
intracavity field anywhere in the detector as a linear combination of basis functions related to the 
eigenmodes of one of the Fabry-Perot arm cavities, and we take great care to preserve numerically 
the nearly ideal longitudinal phase resonance conditions that would otherwise be provided by an 
external servo-locking control system. We have included the effects of nonlinear thermal focusing 
due to power absorption in both the substrates and coatings of the mirrors and beamsplitter, the 
effects of a finite mismatch between the curvatures of the laser wavefront and the mirror surface, 
and the diffraction by the mirror aperture at each instance of reflection and transmission. We 
demonstrate a detailed numerical example of this model using the MATLAB program Melody for 
the initial LIGO detector in the Hermite-Gauss basis, and compare the resulting computations of 
intracavity fields in two special cases with those of a fast Fourier transform field propagation model. 
Additional systematic perturbations (e.g., mirror tilt, thermoelastic surface deformations, and other 
optical imperfections) can be included easily by incorporating the appropriate operators into the 
transfer matrices describing reflection and transmission for the mirrors and beamsplitter. 

PACS numbers: 04.80.Nn, 07.05.Tp, 07.60.-j, 07.60.Ly, 42.60.-v, 42.60.Da, 44.05.-|-e, 95.55.Ym 



I. INTRODUCTION 

In the present decade a number of long-baseline laser 
interferometers are expected to become operational, and 
begin a search for astrophysical sources of gravitational 
radiation. ||, |j These include the Laser Interferom- 
eter Gravitational Wave Observatory (LIGO),|^, ||] the 
VIRGO project, § the GEO 600 project, (t) and the 
TAMA 300 project. H All of these detectors will employ 
a variant of a Michelson interferometer illuminated with 
stabilized laser light. The light will be phase-modulated 
at one or more radio frequencies, producing modulation 
sidebands about the carrier frequency which provide a 
phase reference for sensing small variations of the in- 
terferometer arm lengths. B Gravitational radiation will 
produce a differential length change in the arms of the 
Michelson interferometer, resulting in a signal at the out- 
put port. 

In Fig. |l|, we show the configuration of the initial LIGO 
detector, a power-recycled Fabry-Perot Michelson inter- 
ferometer (PRFPMI). The light from a stabilized laser 
source enters the interferometer, which is comprised of an 
asymmetric Michelson interferometer with Fabry-Perot 
arm cavities. The arm cavities consist of polished in- 
put and end test masses (ITM and ETM, respectively) 
whose coated surfaces also act as mirrors to resonantly 
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FIG. 1: Schematic configuration of the initial LIGO detec- 
tor, a power-recycled Fabry-Perot Michelson interferometer 
(PRFPMI). 



enhance the laser light. An additional power-recycling 
mirror (PRM) placed between the laser and beamsplit- 
ter increases the total light power available to the arms, 
by forming a "recycling cavity" together with the beam- 
splitter and arm cavity input mirrors. The sideband fre- 
quency and interferometer lengths are adjusted so that 
the sidebands resonate in the recycling cavity but not in 
the arms, while the carrier resonates in both sets of cavi- 
ties. This allows extraction of control signals for all four 
length degrees of freedom. 

The noise level in the detection band at the interferom- 
eter output port must be held below the required strain 
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sensitivity. The primary noise sources limiting the in- 
terferometer sensitivity are seismic noise at frequencies 
below 100 Hz, thermal noise between roughly 100 and 
300 Hz, and photon coimting noise at frequencies greater 
than 300 Hz. Photon counting noise is suppressed by 
using sufficient laser power; a design strain sensitivity of 
roughly 10~^^ over a time period of a millisecond will 
require a power of several hundred watts incident on the 
beamsplitter. 

A concern in the operation of an interferometer at such 
high power levels is optical distortion. The presence of 
absorption in the coatings and substrates of the optics 
can give rise to non-uniform heating due to the Gaus- 
sian intensity profile of the laser light. This heating will 
cause a temperature rise in the substrate, which, coupled 
to the temperature dependence of the substrate index of 
refraction, will cause a distortion of the wavefront of the 
light transmitted through the substrate. For example, 
the initial LIGO design provides for fused silica input 
test masses which are 10 cm thick. To ensure the sta- 
bility of the arm cavity resonant mode, the ITM radius 
of curvature is set to be several times the length of the 
arm cavity, or of order 10 km. By comparison, the sub- 
strate heating due to a nominal substrate and coating 
absorption of 5 ppm/cm and 1 ppm respectively causes 
the generation of a thermal lens in the substrate with 
a nominal focal length of 8 km. This significant distor- 
tion can have important consequences for both the power 
buildup of the RF sidebands in the recycling cavity and 
the interferometer sensitivity. 

There have been a number of efforts made to under- 
stand the effects of the propagation of distorted wave- 
fronts through gravitational wave interferometers, in- 
cluding an approximate treatment of wavefront distor- 
tion in power-recycled and dual-recycled systems due to 
imperfections in the shape or alignment of some opti- 
cal components, Jill an approximate treatment of ther- 
mal lensing and surface deformation in early Fabry- 
Perot Michelson designs, and an approximate de- 
scription of thermal lensing in alternate signal extraction 
designs. [|l3[ Q Each of these discussions emphasized the 
effects of optical distortion on a particular subset of the 
optical performance characteristics of an interferometer. 
In principle, fast Fourier transform methods ^ |l^ 
can be used to simulate a thermally loaded interfer- 
ometer, but the computational cost of an iterated self- 
consistent FFT numerical analysis of the corresponding 
intracavity fields would be prohibitive. 

By contrast, we have developed a flexible and robust 
analytical and numerical model of the effects of wavefront 
distortion on the optical performance of a gravitational- 
wave interferometer that requires far fewer computa- 
tional resources than an FFT approach. The model in- 
cludes accurate contributions to the thermal distortions 
in all substrates due to optical absorption in the substrate 
and both coatings of each test mass,[|l8j the effects of 
any mismatch between the curvatures of the laser wave- 
front and the mirror surface at each reflection, and the 



diffraction by an aperture at each instance of reflection 
and transmission. Operators representing these physi- 
cal effects on the reflection and transmission properties 
of mirrors and beamsplitters are explicitly incorporated 
into a set of transfer (or scattering) matrices that relate 
output fields to input fields incident on the surfaces of 
those elements. The corresponding transfer matrices of 
more complex optical components (e.g., a Fabry-Perot 
interferometer) can subsequently be constructed in a hi- 
erarchical fashion from those of simpler systems. Us- 
ing this approach, the model determines a self-consistent 
steady-state solution to the interferometer power buildup 
in the presence of optical distortion, and employs an in- 
ternal numerical mechanism to preserve the nearly ideal 
longitudinal phase resonance conditions that would oth- 
erwise be provided by an external servo-locking control 
system. We have developed computer code representing 
this model that is sufficiently flexible to allow new phys- 
ical effects to be added to the simulations with a mod- 
icum of effort, and that can be stopped, interrogated, 
and restarted at any stage of a simulation. Our imple- 
mentation of the artiflcial resonance-locking mechanism 
provides an exponential reduction in the execution time 
required for a self-consistent simulation by dynamically 
adjusting the frequency response of the interferometer. 

Our current model employs basis function expan- 
sions of the electromagnetic field everywhere in the 
interferometer . p^, p2[ In principle, these ba- 

sis functions can be spatial eigenfunctions of the inter- 
ferometer, but in practice, imperfect mode-matching of 
the multiple cavities comprising the full interferometer 
and the presence of apertures in the system means that 
round-trip spatial propagators are not Hermitian, and 
the eigenfunctions are not power-orthogonal. Q 
Therefore, we choose an expansion of the intracavity field 
in an arbitrary set of power-orthogonal unperturbed spa- 
tial basis functions that are not necessarily eigenfunctions 
of any subsystem of a thermally loaded, perturbed inter- 
ferometer. Since we can move from one spatial represen- 
tation to another simply by recomputing the propagator 
matrix elements, the relative merit of a set of possible ba- 
sis choices is determined by the number of basis modes 
needed to describe the intracavity field accurately. Our 
analysis of the operators encapsulating the physics of the 
perturbative phenomena described in the model is in- 
variant under such a change of spatial basis. Neverthe- 
less, the choice of an incomplete subset of basis functions 
inevitably introduces false optical losses because power 
is transferred to higher-order spatial modes which are 
not included in a necessarily finite numerical simulation. 
Therefore, we introduce unitary approximations for non- 
diffractive operators that monotonically improve as more 
basis functions are included in the simulation. 

In Section || of this paper, we describe our representa- 
tion of the intracavity electromagnetic field, and the op- 
erators incorporated into the transfer matrices of the sub- 
systems comprising components of typical interferomet- 
ric gravitational wave detectors. We develop a unitary 
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approximation for non-difFractive operators expanded in 
finite-dimensional basis sets, and, in our discussions of 
the Fabry-Perot and Michelson interferometers, we intro- 
duce abstract idealizations of the control systems needed 
to maintain different intracavity resonance conditions. In 
Section [II, using these building blocks we construct a 
model of the initial LIGO interferometer, and we point 
out that the basis functions chosen for a matrix descrip- 
tion of the interferometer need not be spatial eigenfunc- 
tions of the fully coupled resonators. As a preliminary 
verification of our method, we compare the predictions 
of our numerical implementation of this model in a set 
of well-defined linear test cases with the output gener- 
ated by a fast Fourier transform computer code suite. 
Finally, we simulate the optical performance of the ini- 
tial LIGO interferometer under full thermal load, using 
a straightforward nonlinear solution method that com- 
putes the steady-state fields everywhere in the system 
after only a few seconds on an ordinary desktop com- 
puter. 



II. MODELS OF OPTICAL COMPONENTS 
AND STRUCTURES 

A. Representation of the Electromagnetic Field 

Throughout this work, we represent a propagating 
laser electric field E(r,t) with angular carrier frequency 
ujQ as the real part of a product of a real time- independent 
polarization unit vector e, a complex amplitude function 
E{r,t), and the carrier wave function exp [i{koz — upt)], 
where Ao = 27r/fco = ^ttc/uiq is the laser wavelength: |23| 

E(r,i) = Rc{eE{r,t)c-Kp[i{koz - LUQt)]} . (1) 

We use the amplitude function E{r, t) to represent any 
closely-spaced frequency components (such as radio- 
frequency-modulated sidebands) as 



E{r,t) =^Eq{r,t)exp [i{kqZ 



,t)]. 



(2) 



where the summation is taken over all frequency compo- 
nents (including the carrier), Eq{y,t) is the complex am- 
plitude function of component 5, and Acjg = Uq — loq = 
Akq/c <^ ujQ is the angular frequency shift of component 
q. By convention, the carrier is labeled by g = 0, and 
therefore Awq = 0. We express all field amplitudes in 
units of ^yWJcm^, so that the time-averaged intensity 
carried by E(r,t) is /(r,t) = ^ \Eq{r,t)\''. 

We will represent the transverse spatial dependence 
of the field amplitude Eq{r,t) at any propagation plane 
within a laser interferometric gravitational- wave detector 
as a linear superposition of a set of amplitude basis func- 
tions with fixed transverse spatial profiles. In general, 
these amplitudes are eigenfunctions of a non-Hermitian 
integral transform equation (generally a composition of 



many consecutive integral transforms, each of which car- 
ries the field from one aperture to the next) that describes 
round-trip propagation through a specific resonant sub- 
system of the unperturbed interferometer. These reso- 
nant eigenfunctions can then be propagated throughout 
the remainder of the interferometer, defining a unique set 
of spatial basis functions. ||2^ Therefore, given the exis- 
tence of a numerically complete set of these transverse 
spatial eigenfunctions, we can expand the sideband am- 
plitude function as 



Eq{r,t) 



(3) 



Hence, we can convert any composition of two consec- 
utive integral transforms into a simple matrix product. 
In practice, we are interested in the steady-state charac- 
teristics of the detector under thermal load, and we will 
clearly indicate the propagation plane under discussion, 
so we will usually suppress the explicit coordinates {z^t] 
in later sections. 

As an example, consider the simple case of laser field 
propagation through a distance L = Z2 — zi in vacuum. 
Suppose that we choose to represent the spatial and tem- 
poral properties of this field using N transverse eigen- 
functions Wm„(r) and Q frequency components, includ- 
ing the carrier and both positive and negative sidebands 
for each nonzero radio modulation frequency. We arrange 
the expansion coefficients into an NxQ ordered matrix E 
with elements Emnq, where each row represents the spa- 
tial eigenfunction corresponding to a particular choice of 
{to, n} and each column corresponds to a particular fre- 
quency component q, and we assume that L has been 
chosen to provide resonant excitation of the fundamental 
carrier spatial mode Eqqq. In this basis, the spatial con- 
tribution to the vacuum propagation kernel is given by a 
N X N diagonal matrix G with elements 



G 



mn.m'n' 



(4) 



where A ipmn = Vmn — is the net Gouy phase (rel- 
ative to the fundamental spatial mode) accumulated by 
mode mn over the propagation distance L, and the net 
longitudinal optical path length contribution (relative to 
the carrier) is given by a Q x Q diagonal matrix U, with 
elements 



Vtqqi = exp(iAwqr) 5, 



qq' 1 



(5) 



where t = L/c\s the propagation time. The components 
of the field following the propagation step can then be 
computed using the sparse matrix product 



E{z2) = GE{zi)n. 



(6) 



We are particularly interested in the more complex ex- 
ample of a power-recycled Fabry-Perot Michelson laser 
gravitational-wave interferometer. In its unperturbed 
state (i.e., in the absence of aperture diffraction and 
thermal wavefront distortions), this detector configura- 
tion is described well by Hermite-Gauss basis functions 
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under perfect mode-matching conditions. The infinite- 
aperture Hermite-Gauss basis is orthonormal (since it 
is generated by a Hermitian propagation equation), so 
that field amplitude expansions using this basis are 
power orthogonal, allowing time-averaged intensities 
to be calculated using the sum of the inner products 
\Eg{r,t)\^ = -^mn \Emnq{z,t)\\ If both Gouy and op- 
tical path length phase contributions are accumulated 
using Eq. (^, then we can represent an arbitrary non- 
astigmatic Hermite-Gauss basis function Umn{x,y) at a 
given reference plane using 



i{x,y) = 




X iJ, 



X exp I— ^ 



(7) 



where Hn{x) is the Hermite polynomial of order n, w 
is the beam spot radius, and R is the wavefront radius 
of curvature at that reference plane. Throughout this 
work, we will define the x-z plane to be "horizontal," 
containing the central beam paths of the two Fabry-Perot 
arm cavities in Fig. |^, and the y-z plane to be "vertical," 
perpendicular to the horizontal plane and parallel to the 
local gravity gradient of the Earth. 



B. Mirrors 

1. Aperture Diffraction and Mirror-Field Curvature 
Mtsmatch Operators 

Consider the general problem of the reflection of a 
propagating electromagnetic field by a perfectly reflect- 
ing cylindrically symmetric mirror M with spherical ra- 
dius of curvature Rm and aperture diameter 2a. The 
reflection can be described in a single step taken from 
the aperture A at the reference plane z^ conventionally 
located just before the reflection from A4 to the field po- 
sition z'' just after reflection has occurred. If the trans- 
verse field at is E{x' ,y' , z^), then the corresponding 
field at z'' is given by Huygens's integral in the Fresnel 
approximation asp^ 

E{x,y,z>)= f dx'dy' K{x,y;x',y')E{x',y',z<), (8) 

J A 

where the functional form of the forward propagation 
kernel for a cylindrically symmetric mirror is |22| 



K{x,y;x',y' 



exp 



TT 2 



Ao Rm 
X S{x - x')6{y - y'). 



2 , ,2 

y 



(9) 



Here the leading negative sign accounts for the Maxwell 
boundary condition at the perfectly reflecting mirror sur- 
face. 



If we expand the field at both reference planes using a 
numerically complete set of cartesian basis functions as 
m Eq. (|), then we can compute the matrix elements of 
the forward propagation kernel in that basis as 



mn^m'n' 



dxdy II dx' dy' K{x,y;x' ,y') 



A 



X ul^n{x,y,Z^)u„i'n'{x',y',Z<) (IQ) 



dx dy exp 



An R 







(x' + y') 



X ''^Inni-^T Z^) ''^m'n' {x, y, z'^) 

Here the function u]^^{x, y, z^) is not generally the con- 
jugate transpose of Umn{x,y, z^). Instead, these basis 
functions represent fields propagating through the system 
in the reverse direction, and are obtained by solving the 
Huygens-Fresnel integral Eq. (^) using K^{x, y; x' , y'), 
the transpose of the forward-propagation kernel. |2^, 

Although the refiection-diffraction problem is specified 
completely by Eq. (|l^), it can introduce inconsistencies 
in practical implementations of optical resonator simula- 
tions. For example, when a finite number of basis func- 
tions are used to represent the propagating field, the ma- 
trix calculated using Eq. (10) generates two sources of op- 
tical loss: genuine diffraction loss due to field truncation 
by the aperture, and power transferred to higher-order 
optical modes which are ignored in the simulation. The 
latter contribution is unphysical: in the limit a — > oo, 
when a numerically complete set of basis functions is used 
the propagation matrix must conserve energy (i.e., must 
be unitary.) Therefore, for the purposes of our simula- 
tions we will explicitly allow diffraction loss due to aper- 
ture truncation only, and separate the effects of aperture 
diffraction and reflection into two consecutive propaga- 
tion steps using the formulation 



K^CA, 



(11) 



where A is a matrix operator representing purely trans- 
missive aperture diffraction, and C is a unitary matrix 
operator describing direct focusing by a spherical mirror 
at normal incidence. 

As a concrete example, we first calculate the matrix 
elements of A in the Hermite-Gauss basis, treating the 
truncation as a perturbation and describing the field both 
before and after the aperture using a linear superposition 
of a flnite set of unperturbed Hermite-Gauss functions. 
Therefore, uj„„(x,?/) = uj„„(a;,y), and we define the an- 
alytical integral 



Amn.m'n' ~ / / dx dy U^j^{x , y) Um'n' {x , y) 

JJa (12) 

where a = 2(a/w;)^, and the integration is performed 
over the circular aperture defined by \/ x"^ + y'^ < a. We 
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specify a basis set by designating the transverse modes for 
a given frequency component q in order of increasing m + 
n first, and then in order of increasing n. For example, if 



we use a Hermite-Gauss basis with m+n < 2, then for the 
column vector Eg = {E'oog, £^iog, £^oig, E20q, Eu^, Eo2q}'^ 
we obtain the matrix 
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(13) 



Note that nonzero elements of I{a) satisfy the relation 
/™„,™'n'(a) ~ 0[a^('"+"+'"'+"')]. 

Next, we ignore aperture diffraction and calculate the 
matrix elements of C in the Hermite-Gauss basis. Strictly 
speaking, the matrix elements of C can be calculated 
using Eq. (10) in the limit a — > cxd. However, as described 
above, we seek an explicitly unitary approximation for 
C that monotonically improves as we include more basis 
functions in the field expansion. When Eq. (|l^) is applied 
to a Hermite-Gauss basis function with field curvature 
parameter Rp, this parameter is transformed according 
to the ABCD propagation law as ^/R'p = 1/Rf — 2/Rm- 
In the unperturbed case where Rm ~ Rf, the curvature 
parameter of the backward-propagating basis function is 
therefore —(—Rp) ~ Rp, and using Eq. (|^) we have 

C = exp(z7c), (14) 



where 



and 



nw" I 1 

Rp 



A 



1 

Rm 



dx dy [x^ + y^) 



(15) 



(16) 



y)\- 



As the number of Hermite-Gauss basis functions used to 
represent Eg grows, exp(i7c) C, and since c is sym- 
metric, exp(i7c) is explicitly unitary. Using the gener- 
ating function of the Hcrmitc polynomials, we obtain the 
matrix elements 



2 

n,n' ^ 



(17) 



where 



+ - V'(m' + l)(m' + 2) 5ra,m'+2 



(18) 
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FIG. 2: The fractional error in [A 'nn.nnl obtained by comparing 
the exact value given by Eq. ([L9|) with the TEMoo element of 
the approximate propagator K = C A, where A is provided 
by Eq. (|l|) and C by Eq. (111. 



As a straightforward test of our unitary approxima- 
tion, we can compare the separated propagator given by 
Eq. ( |rT| ) with the exact propagator specified by Eq. ( |lQ ) 
for the TEMqo Hermite-Gauss mode. An explicit calcu- 
lation using Eq. (Ilfl) yields 



1 _ g-(l + j7)c 



00,00 — 



«7 



(19) 



-I- -^J{m + l)(m -I- 2) 6m',m+2- 



In Fig. ^ we plot the fractional error in li^oo.ool accrued 
by computing K = C A, where A is provided by Eq. ( p^ ) 
and C by Eq. (|^. For an incident field with a spot 
radius w = 3.64 cm and a wave front radius of curva- 
ture Rp — 10 km at a wavelength A = 1.0642 /im, we 
show relative errors for different values of the mirror's 
cylindrical substrate radius and radius of curvature as a 
function of the maximum mode index N . The number 
of basis modes with m + n < N used in the calculations 
of A and C is given by ^{N + 1){N + 2). As we expect, 
the relative error decreases significantly as the number of 
basis modes increases. 
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2. Thermal Focusing Operators 

In an interferometric gravitational-wave detector, the 
test masses (as well as any optics which form power 
and/or signal recycling cavities) are thick cylindrical 
blocks of a transparent material, such as pure fused sil- 
ica. A high-reflectance coating (HR) is applied to the 
broad inner face of the cylinder, forming a mirror; if 
transmission through the mirror substrate is necessary, 
then an antireflecting coating (AR) will be applied to the 
outer surface as well. Under test, then, laser power may 
be absorbed in the coatings and/or the mirror substrate 
(SS), and the corresponding heat flow into the substrate 
will result in an inhomogeneous temperature distribution 
throughout the mirror. This temperature gradient will 
cause a refractive index gradient to form, converting the 
substrate into a thermally-generated thick lens. 

In general, the propagation of the laser field through 
the substrate will be described by an integral trans- 
form equation that incorporates the effects of position- 
dependent optical phase shifts caused by the thermal 
load. As discussed above, by choosing a suitable set of 
unperturbed basis functions (i.e., in the absence of ther- 
mal and mechanical perturbations) we can convert the 
integral transform equation into a scattering matrix op- 
erator which redistributes energy contained in an initial 
configuration of basis functions into the final set of cor- 
responding functions that emerges from the substrate. 
In order to determine the elements of this matrix oper- 
ator, we must first find the temperature distribution in 
the substrate due to absorption in both the coatings and 
the substrate, and then we must compute the matrix el- 
ements of the longitudinal optical path-difference (OPD) 
phase shift introduced by the temperature gradient. 

Hello and Vinet||l^ have calculated both the steady- 
state and transient temperature distribution throughout 
the simplified mirror shown in Fig. ^ for the case where 
the change in the temperature due to the absorption of 
optical power remains small compared to the external 
ambient temperature Tq. This assumption linearizes the 
radiative boundary conditions at the surfaces of the mir- 
ror, and allows the temperature increase due to coating 
and substrate absorption to be calculated independently 
(and then summed) for some incident power P. As shown 
in Fig. ^, they treated the mirror substrate as a thick 
disk spanning the cylindrical coordinates < r < a and 
— h/2 < z < h/2, and the coating is concentrated in an 
infinitesimally thin layer a,t z — —h/2. After an alge- 
braic simplification of Hello and Vinet's steady-state re- 
sults, we find that the temperature distributions through- 
out the substrate due to absorption in the substrate and 
coating are, respectively, 



Ts{r,z) = P 



UsCL 



E 



1 - 2TAk cosh Cfc - 
a 



(20) 



/ \V,Coating\ 




Substrate 



z = -h/2 z = +h/2 

FIG. 3: Schematic diagram of a singly-coated mirror and the 
corresponding coordinate system used by Hello and Vinetjist. 
Both the mirror and the incident power distribution are as- 
sumed to be azimuthally symmetric. The total incident power 
P can be absorbed in the coating and /or the substrate, re- 
sulting in a total temperature distribution T{r, z) and a cor- 
responding thermal lens in the substrate. 



Tc(r,z) = P'^^Pk 



Ak cosh 



Bk sinh C 



(21) 



where is the power absorption of the coating, ag is the 
absorption coefficient of the substrate material. 



Bk = 



1 



2 [Cfe sinh (7fe) + T 
1 



cosh (7fc)] ■ 



2 [Cfe cosh (7fe) -f r sinh (7^)] ' 



(22a) 
(22b) 



7fe = Cfeft,/2a, T = AeaTQa/kT, tr is the Stefan-Boltzmann 
constant, e is the total spherical emissivity of the sub- 
strate, and kx is the thermal conductivity of the sub- 
strate. For a — 12.5 cm (consistent with the value cho- 
sen for the LIGO test masses) and Tq = 300 K, we have 
T = 0.277. The terms in these series are characterized by 
the roots of the equation 



CJi (C) - rJo (C) = 0, 



(23) 



which are reasonably well-approximated by the zeros of 
Ji (C), Cfc = (fc + l/4)7r for k € {0, 1,2,.. .}, when r < 1. 

The constants pk are the coefficients of a Dini series 
expansion jlSf of the azimuthally symmetric incident in- 
tensity distribution /(r), given by 



Iir) = Y,PkJo{Ck^) 

k 



(24) 



A straightforward calculation using Eq. (|2^) yields 



Pk 



{ek + r')Ji{<:k) 
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TABLE I: Material constants for fused silica used in our sim- 
ulations. 



Constant 



Name 



Value 



Units 



rj Refractive index 

drj/dT First-order rj-T coefficient 

fcr Conductivity 

e Total spherical emissivity 

Us Bulk absorption coefficient 

(Ts Bulk scattering coefficient 0.5 x fO 



i. 44968 

8.7 X fO"'' K"^ 

f .38 W/m K 
0.76 

5 X iO"'' m-i 



If we assume that the field is primarily comprised of 
the fundamental {00} mode, then we can neglect the 
slight heating arising from higher-order modes. When 
the incident intensity describes a unit power TEMqo 
Gaussian beam with spot radius w, we have /(r) — 
(2/7rw^) exp (— 2r^/t(;^) . Furthermore, if w/a <C 1, we 
can extend the upper limit of integration in Eq. ( p5| ) to 
infinity, thereby obtaining the analytic expression 



Pk 



exp(-/3fc). 



(26) 



where Pk = |(Cfcw/a)^ 

The z dependence of the substrate temperature dis- 
tribution given by Eq. (^0|) is typically very weak when 
w/a <C 1 and the boundary conditions are radiative. 
In this case, near r = we can use the identity J^kPk — 
/(O) = 2/7rw^ to obtain the approximate temperature 
distribution Tsir) - Ts{0) ^ {apPs/2'KkT){r /wf . In the 
thin lens approximation, this quadratic temperature gra- 
dient yields a weak focal length 



/ 



(27) 



ashP dri/dT 

The material constants for fused silica are shown in Ta- 
ble |. If we assume typical LIGO parameter values for 
the test masses, then we have w = 4 cm, h = 10 cm, and 
P = 500 W, giving / = 16 km, a length only four times 
greater than that of a LIGO Fabry-Perot arm cavity. 

The thermal load on the substrate causes a position- 
dependent phase shift across any wavefront propagating 
through the mirror. Defining the longitudinal optical 
path-difference (OPD) phase (f> as 

cj^ir)^--^ dzT{r,z), (28) 

we obtain for the substrate and coating contributions 
27ra^ drj sr-^ pk 



Mr) = Ps 



Mr) = Pc 



XakTdT^Cl 



2TAk 

Ik 



sinh (7fe) 



/o (a 



(29a) 



27ra^ dr\ -^—^ pk 

x^df^a 



2 Ah sinh 



(7fe) Jo {Ck 9 > (29b) 




HR 



A 
Z 



AR 



FIG. 4: Schematic diagram of field distributions in a mirror 
of thickness h, consisting of a substrate (SS) with applied 
high-reflectance (HR) and antireflecting (AR) coatings. 



where Ps = aphP and Pc = UcP are the powers absorbed 
in the substrate and coating, respectively. If the thick- 
ness of the substrate is small compared to the wavefront 
radius of curvature, then we can approximate the propa- 
gator describing transmission through the substrate as a 
diffraction step through the aperture followed by a direct 
integration of the OPD across the aperture to describe 
the accumulated distortion. If we define the "edge-to- 
center" OPD phase as A(j>{r) = (t){r) — ^(0), then we can 
compute the total substrate thermal distortion matrix 
operator using 



Smn,m'n' = J J dx dy cxp [iA</)(r)] 

X 'uLn{x,y)Um'n'{x,y) 



(30) 



for both substrate and coating absorption. 

As wc discussed above for the wavefront-mirror cur- 
vature matrix C given by Eq. (p^, we seek a unitary 
approximation of the substrate thermal distortion ma- 
trix to ensure that power is conserved even when a finite 
set of unperturbed basis functions is used to compute 
the matrix elements. Once again, we begin by comput- 
ing the matrix elements of the exponential arguments 
A(/)s(r) and Acpdr) as 



s-.mn.m' n' 



Ps 



Eh. 



2na^ drj 



XokrdT^Ck 



Ik 



■ 2Ak sinh (7^ 



mn.m' n 



<31a) 



c-.mn.m' n' 



Pr 



Pk 



2710? df] 



Xokr dT ^ C,k 

k 



E 



X 2AfeSinh(7fc)e-'^''-/fc; 



mn.m' n' 



(31b) 
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where we have defined the matrix integral Ik as 

oo 

4;mn,m'«' = e'^'' J J dx dy Jq (Ck ^) " 

— OO 



(32) 



Since the integrand is zero at r = + = 0, our def- 



inition of A(/>(r) allows us to explicitly ignore any phase 
shift accumulated due to the peak thermal change in 
the substrate refractive index, and instead we include 
only the contribution to the wavefront distortion arising 
from the shape of the temperature distribution. Using 
the same example Hermite-Gauss basis described in Sec- 
for m + n < 2 and a given value of fc, we 



tIBl 



tion 

obtain the square matrix elements 



IkilSk) 



1 




1 



71^^ 



Pk 


























4Pk 



Jf-'k 








2/3fc- 







4Pk 



1 - 2/3, + I/32 



(33) 



Note that nonzero elements of Ik{(3k) satisfy the relation 

Collecting results, we can compute the net effective 
substrate OPD matrix operator in our unitary approxi- 
mation as 



5 = exp{i [P,$, + (P,, + P,)$,]}, 



(34) 



where Pg, Ph, and Pa are the total TEMoo optical powers 
absorbed in the substrate, HR coating, and AR coating, 
respectively. The power distribution throughout the mir- 
ror volume is shown schematically in Fig. ^ for a dual- 
coated cylindrical mirror with substrate thickness h. A 
laser field is incident on the HR coating (the "front" 
of the mirror), and another field is incident on the 
AR coating (the "back" of the mirror). Since we have 
assumed that the field is primarily comprised of the fun- 
damental {00} mode, then the total power absorbed in 
the coatings and substrate depend only on the field am- 
plitude expansion coefficients -Foog i as listed in Table ||. 
Here ahr and aar are the dimensionless fractional power 
absorptions of the HR and AR, respectively, and is the 
power absorption coefficient of the substrate, with units 
of inverse length. The reflected field E~ is the superpo- 
sition of the prompt reflection of F+ and the residual of 
F" which has been transmitted through the HR. 



3. Transfer Matrix 

We can extend our analysis of the reflection and trans- 
mission processes and build a general transfer matrix for 
fields incident on the mirror from either direction. First, 
we choose the laser-industry standard conventions for the 
phases of the amplitude reflection and transmission co- 
efficients for a quarter-wave dielectric stack. Therefore, 
the Maxwell boundary conditions for the HR generate 



TABLE II: Total absorbed TEMoo powers for each of the 
three regions shown in Fig. ^ summed over the lowest-order 
transverse modes of all active sidebands. The factors of i 
arise from our normalization convention for the laser electric 
field. 



Region 




Absorbed Power 




SS 
HR 
AR 


Ps 

Ph 

Pa 


= «./^5E,(lSo'o,r + 
= E, (|i^o+oJ' + 
= «-^E, (1^^0+0,1' + 


\Foof) 
Foof) 



the amplitude coefficients —r for reflection and i t for 
transmission, regardless of propagation direction. Losses 
in the HR due to absorption and scattering are included 
implicitly when r"^ + < 1. However, we must explicitly 
account for losses due to transmission through the sub- 
strate and AR coating using the amplitude coefficient 



s)h/2 



\/\ - {aar + Sar), 



(35) 



where ag is the scattering loss per unit length in the 
substrate, and Sar is the scattering loss in the AR coating. 

The propagator for reflection from the vacuum side of 
the HR is given by Eq. (^), and the resulting perturbation 
matrix for a mismatch between the radius of curvature 
of the mirror and the wavefront is given by Eq. ( p^ as 
(7(7) = C~ . Since the propagator for reflection from the 
substrate side of the HR is also given by Eq. (|) with 
—2/Rm 2n/i?M,[l3) where n is the refractive index of 
the substrate, the corresponding perturbation matrix is 
C{—n^) = . Similarly, the propagator for transmis- 
sion through the HR from the substrate to the vacuum 
is given by Eq. (|l|) with -2/Rm ^ {n - 1)/Rm,^ so 
the transmission curvature mismatch matrix operator is 
CKn- 1)7/2] = (C-C+)V^ 



In Section piD2 and Section IIIB, we will adjust the 
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FIG. 5; Schematic diagram of the generalized mirror transfer 
matrix given by Eq. (pq) . The HR-coated surface of the mirror 
is on the left. The labels and 2^ denote the input and 
output reference planes, respectively, near the HR surface, 
while (^^ and are the input and output reference planes 
near the AR surface. 



microscopic positions of some of the mirrors to ensure 
that the pertmbed interferometers maintain the appro- 
priate resonant phase conditions as the system evolves. 
If the mirror in Fig. ^ is moved in the positive z direction 
(i.e., toward the right in the diagram) by the microdis- 
placement tS.z < A, then reflection from the "front" of the 
mirror introduces a relative phase exp(-|-i2A;A2;), while 
reflection from the "back" of the mirror requires a phase 
adjustment of exp(— i2fcAz). 

Collecting these results, the flnal mirror transfer ma- 
trix depicted in Fig. ^ is 



'e-' 




't- 


R-' 






E+ 




R+ 


T+ 




F+ 



(36) 



where the transmission operators are 



T- = itt,{C-C+Y''^ S A, and 



r+ = ittsS{c+c-f''^ A 



and the reflection operators are 

-re+^^fe'^^C^ A, and 



R- 
R+ 



(37a) 
(37b) 



(38a) 
(38b) 



In the lossless case where + — t1 — 1 and A is 
the identity matrix, the explicit unitarity of the , C~ , 
and S matrix operators guarantees that + \ E^\'^ ~ 



C. Beamsplitter 

1. Aperture Diffraction and Thermal Focusing Operators 

The analysis of the aperture diffraction and thermal 
focusing operators for the case of a 45° beamsplitter is 
similar to that of a mirror. However, even though the 
beamsplitter has the same cylindrically symmetric sub- 
strate as the mirror shown in Fig. ^, the non-normal angle 
of incidence will require the elements of the perturbative 
matrix operators to be computed numerically. We will 
assume here that the HR-coated surface of the beamsplit- 
ter is flat, allowing us to ignore both curvature mismatch 
and astigmatism when we build the transfer matrix in 



Section II C 2 



The aperture diffraction operator for reflection from 
the HR-coated surface of the beamsplitter is given by 



A. 



v) {x,y), (39) 



where £ represents an elliptical aperture with semimajor 
(y) axis a and semiminor (x) axis a/^/2. Strictly speak- 
ing, the diffraction operator for transmission through the 
substrate is more complicated than Eq. ( |39| ) because of 
refraction at each vacuum-dielectric interface, particu- 
larly since the effective clear aperture (in the horizon- 
tal plane) after propagation through the substrate has 
been reduced by the distance h/2ri, or about 1.4 cm for 
LIGO. In principle, we can construct separate operators 
for reflective and transmissive aperture diffraction, but 
in practice we have found that the resulting effect on the 
recycled power enhancement of the standard PRFPMI 
configuration shown in Fig. |l| was negligible. Hence, in 
our simulations, we have made the simplifying approxi- 
mation that Eq. ( ^9| ) can be used to represent aperture 
diffraction effects for both reflection and transmission. 

The distribution of absorbed power within the beam- 
splitter volume is illustrated in Fig. ^, where we have in- 
troduced the primary propagation paths labeled parallel 
(II) and perpendicular (_L) relative to the propagation of 
the input laser fi eld in Fig. ^ Following the conventions 
of Section II B 2 as closely as possible, we then label the 
input fields incident on the front (HR-coated) side of the 
mirror with sign(k • z) > as F+H and F^-^, and those 
incident on the back (AR-coated) side of the mirror, with 
sign(k • z) < as F"" and F~'^. Together with the four 
corresponding output fields shown in Fig. ^, these fields 
generate five distinct optical power absorption regions. 
The total TEMqo optical power absorbed in each of these 
regions is listed in Table III, where h' = h/ ^Jl — \/2rf 



is the distance traveled through the substrate by the re- 
fracted beams. 

Diffusion of heat from each of the optical absorption 
regions contributes to the total OPD distortion for both 
the parallel and perpendicular propagation paths. Fo l- 



lowing the conventions we established in Section II B 2 
and Eqs. (pT|), for the parallel propagation path shown 
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can be written using our unitary approximation as 



FIG. 6: Interaction of the two primary propagation paths — 
parallel (||) and perpendicular (_L) — with the five optical 
absorption regions that cause thermal focusing in the beam- 
splitter substrate — SS", SS"^, HR, AR", and AR-^. 



exp< I 



p! $f + p„-^$[! 



exp< I 



-Ph K + Pi < 



P^ * 



and (40a) 



(40b) 



For the purpos es of the initial LIGO simulations dis- 
cussed in Section III E , we first used finite-element soft- 
ware to compute the temperature distribution every- 
where in the beamsplitter substrate due to absorption 
of 1 Watt of TEMqo optical power (choosing the initial 
LIGO spot size at the beamsplitter) in each of the five re- 
gions shown in Fig. |6|. Fo llowing the same basic approach 



as that of Section II B 2 , we then computed the matrix 
elements of the five parallel OPD operators in Eq. (|40| ) 
using a Hermite-Gauss basis having the same spot size. 
In our simulation runs, we simply scaled each matrix by 
the appropriate absorbed power, and then computed S"" 
and by matrix exponentiation. 



TABLE III: Total absorbed TEMm powers for each of the five 
absorption regions shown in Fig. y, summed over the lowest- 
order transverse modes of all active sidebands. The factors of 
i arise from our normalization convention for the laser electric 
field. 



2. Transfer Matrix 



Region 



Absorbed Power 



ss" 
ss-^ 

HR 

ARll 

AR^ 



nil 





2 






+ 






2 




( -^(wt 


-f 


P00q\ ) 



F+\\\ 

00<j| 


2 1 
+ 


OOq 


^QOq 1 


2 1 
-1- 


Eoot\ j 


\^OQq 


|2 


\EooI\ j 




1 2 

+ 


rooq 1 , 



in Fig. 1^, we denote the substrate OPD phase operators 
corresponding to the five absorption regions SS", SS^, 
HR, ARII, and AR-^ as $!, $1, and respec- 

tively. Since the perpendicular and parallel propagation 
paths are symmetric under reflection in the x-z plane 
about the z axis, the same set of operators can be used 
to represent the substrate OPD distortion for perpendic- 
ular propagation, but in the order <I>^, $|, <I>^, and 

<I>c. Hence, in a particular basis, the net substrate ther- 
mal OPD matrix operators for the two propagation paths 



As in the case of the mirror described in Section [I B 3 , 
we can construct a generalized transfer matrix for the 
beamsplitter that includes the effects of aperturing, ther- 
mal focusing, and longitudinal microdisplacement. Once 
again, we adopt the industry-standard convention for 
the phases of the quarter-wave amplitude reflection and 
transmission coefficients of the HR coating, and we ap- 
ply the substrate/ AR amplitude transmission coefficient 
given by Eq. (35) with h h' . However, as shown in 
Section [I E 2| , the definition of the beamsplitter position 
is more subtle than that of the mirror: it consists of a 
static contribution z that represents the location cho- 
sen for the beamsplitter when used in cold, unperturbed 
systems, and a dynamic contribution Az that is used 
in heated interferometers to maintain a given intracav- 
ity phase condition. Hence, if the beamsplitter in Fig. ^ 
is displaced in the positive z direction (i.e., away from 
the HR coating along the beamsplitter axis of cylindrical 
symmetry) by the total displacement z + Az, then re- 
fiection from the "front" of the beamsplitter introduces 
an additional relative phase exp[+i2^/2k{z + Az)], while 
reflection from the "back" of the beamsplitter requires a 
phase adjustment exp[—i2^/2k{z + Az)]. 

Collecting these results, the final beamsplitter transfer 
matrix depicted in Fig. |^ is 



, (41) 



"p-ll" 




"T-II 
































E+W 










T+ll 


R+W 




F+W 
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where, in the absence of curvature mismatch, the trans- 
mission operators are 



T^ll = r+" =itt',S^^A, and 
and the reflection operators are 



(42a) 
(42b) 



(43a) 



= -r^;2e-'^^fc(^+^^)s•ll S'-LA, and (43b) 
= -rt;2e-»v^fc(^+^^)s'-L5" A (43c) 



In the lossless case where 



= <'2 = 1 and A = 1, 



the explicit unitarity of the S" and S*^ matrix operators 
guarantees that \E+-^\'^ + iS^Hp = \F+-^\'^ + iF^Hp and 

|£;+ii|2 + |£;--L|2 = |F+ii|2 + |F--L|2. 



D. Fabry-Perot Interferometers 

The propagation schematic diagram of a Fabry-Perot 
interferometer (FPI) is shown in Fig. ^. The front (HR) 
input and output ports of mirror A^i and are coupled 
through a propagation distance L13. We can simplify 
later calculations of the optical performance of a power- 
recycled gravitational-wave interferometer if we first de- 
termine a transfer matrix for the FPI, and then establish 
a procedure maintaining an optimum level of intracavity 
field amplitude enhancement in the presence of system- 
atic perturbations. Although we develop the FPI transfer 
matrix in a basis-independent fashion, and postpone the 
discussion of the choice of an unperturbed basis set until 
we calculate steady-state fields in Section [II C , we note 



that we will use the unperturbed Hermite-Gauss eigen- 
modes of one of the Fabry-Perot interferometers in Fig. |l| 
to develop a basis set for the entire gravitational-wave 
detector. 



1. Transfer Matrix 

Consider the schematic representation of the enhance- 
ment, reflection, and transmission operators of the Fabry- 
Perot interferometer shown in Fig. If we express these 
operators using the transf er mat rices of mirrors M 1 and 
A^3, described in Section [IB3, then the self-consistent 



forward propagation steps between the plane of incidence 
of All at zf- and the plane of incidence of M3 at can 
be represented as 



3q 



^-Gi3 (i?3- F+ + T3- F-) , and (44a) 

(44b) 



^G3i {R^F^^+T^F, 



+_L 



'-II 




(a) Cross-coupling for -|- 1. and 



'-_L 




(b) Cross-coupling for 



and 



FIG. 7: Schematic diagram of the generalized beamsplitter 
transfer matrix given by Eq. (^). The HR surface of the 
mirror is on the upper left. 



where F^ is the amplitude of the external sideband fleld 
q incident on A4i, F^ is the amplitude of the external 
sideband fleld q incident on A^3, F^^ is the amplitude of 
frequency component q incident on reference plane z^, representing these operators are identical.) If the length 



and i^3^ is the corresponding fleld incident on reference 
plane z^ . Here G31 is the Gouy operator describing the 
forward propagation step from zj* to z^, while G13 is 
the corresponding operator describing propagation from 
z^ to . (In a power-orthonormal basis, the matrices 
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1 M 



> 
< 



< 



Z3 



> 



Z3 
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FIG. 8: Schematic diagram of the enhancement, reflection, and transmission operators of a Fabry-Perot interferometer, defined 
by Eq. @, Eq. (||), and Eq. (H), respectively. 



of the vacuum separating A4i from AI3 is L31, then the 
single-pass propagation time is T31 = L^i/c = T13 
Using the concatenated column vectors [E^^ 
and [Fg' F^]'^ , we can construct a system of sparse ma- 
trix equations which can be efficiently solved numerically 
using Gaussian elimination: 



Alternatively, for a purely analytic calculation, we can 
define the global enhancement operator 



H{Auj)^ 1-G{Alu)R- G(Aw)r-, (48) 



1 - GiAujg) R' 
where 



T~ = 









G{Alu,)T' 



and R ~ 



i?3 
R^ 



(45) 



(46) 



and then solve Eq. ( |45|) analytically to obtain 













H^{Auq 



(49) 



and 



G{Alo) = 



e'^'^^'^^G-ii 



where the operators representing the enhancements at 
(47) each of the reference planes in Fig. @ are given by 



i7f [Auj) 
H+{Auj) 



(1 - e''^""^^Gi3 i?3" G31 Ri) e'^'^^"Gi3T3-, 



(1 



„i 2 Aw ri; 



Gi3 i?3 G31 Ri 



J 1 Aw Ti; 



Gi3 i?3 G31 , 



ii^{Au:) = (l-e*2Aw.3iG3ii?-Gi3i?3-) 'e*2^"""G3ii?rGi3r3^, and 

1 



H+iAco) 



2 Aw Tai 



G31 i?]^ G13 i?3 



,i Aw rsi 



G31 Ti 



(50a) 
(50b) 
(50c) 
(50d) 



From Fig. ^, we see that each of the output fields 
and E'^ can be expressed as the sum of a prompt re- 
flection of the corresponding input field and an enhanced 
intracavity field transmitted through the adjacent mirror, 
giving 



E- = T+F^+R+F+ and 



E^„ 



^3 ^3q 



or, from Eq. (|48|) and Eq 

R+ +T+ H{AiUg) 



(51a) 
(51b) 



(52) 



where 

f+ = 



T+ 
T.+ 



and i?+ = 



R+ 
Rt 



(53) 



We can now construct an explicit transfer matrix similar 
to Eq. ( p6[ ) for the FPI as a monolithic optical element. 
Using Fig. || as a guide, Eq. (52) gives: 



TppiiAuJg) Rppj^{Aujg) 

R+P,{AiUg) Tp+pj(A^,) 



(54) 



where the transmission operators are 

T^^AAlu) = T+77f(Aw), and (55a) 
T+H+ {Acu), (55b) 



Lppi(Aa;) 
T+p,{Au;) 
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and the reflection operators are 

i?ppj(Aw) = R+ +T+ H+{Auj), and (56a) 



i?+Pi(A^) = R++T+H-iAL,). 



(56b) 



2. Idealized Simulation of Servo-Controlled Resonator 
Length Locking 

Suppose that the FPI has been initiahzed using a par- 
ticular set of configuration parameters for mirrors A4i 
and M3, the initial thermal loads are negligible, and the 
initial microdisplacements of both mirrors are set to zero. 
At some later time, a perturbation is introduced, such as 
a different curvature for Mi, that alters the resonance 
condition of the cavity. We assume that a biorthonor- 
mal set of unperturbed basis functions itmn(r) has been 
chosen to represent both the intracavity and extracav- 
ity transverse laser fields prior to any perturbations, and 
we wish to compute the corresponding perturbed trans- 
fer mirror matrix given by Eq. (|3^) in that unperturbed 
basis. 

In any experiment using an FPI, the fields emerging 
from the resonator will be monitored and the cavity 
length will be adjusted to maintain optimum enhance- 
ment. However, rather than simulate a servo- mechanical 
control loop, we choose to compute a microdisplacement 
Azi for Ml directly from a numerical determination of 
the net round-trip phase accumulated by the lowest-loss 
carrier eigenmode of the resonator. Beginning and end- 
ing at the reference plane z^, the perturbed round- trip 
propagator at the carrier frequency for the FPI shown in 
Fig. H is 



Kfpi{Azi) 



-i2kAz 



^Gi3 i?3 G31 Ri 



(57) 



where we have chosen to extract from Eq. (38a) the ex- 
plicit phase offset arising from a nonzero value of Azi . 

We can then compute the spectrum of eigenvalues 
of Kppi{0), and select the eigenvalue Aq = lAole"^" 
with the largest magnitude (i.e., the smallest round-trip 
loss). Clearly, the corresponding eigenvector Eq is also 
an eigenvector of the matrix operator Kppi{Azi), with 
the eigenvalue Ai = Aoe~*^'^'^^i. Hence, if we choose 



Azi — — , 
2k' 



(58) 



then the net round-trip phase accumulated by Eq will be 
zero, and that mode will satisfy the resonance condition. 
In our simulations, we maintain a record of the most 
recently chosen lowest-loss eigenvector, and in the case of 
degeneracy we choose the eigenmode having the largest 
value of the inner product with the previous eigenvector. 

Suppose that we have followed this procedure and 
computed the eigenvector Eq of the matrix operator 
A'fpi(0o/2 fc), representing a maximally resonant field 
amplitude at reference plane in Fig. pi We can mode- 
match this field by choosing input fields that optimally 



couple to Eq through the global enhancement opera- 
tor H{0). For example, if we solve the matrix equa- 
tion Eq = H+{0)F+ for F^, where H:^{0) is given by 
Eq. ( ^Ob| ), we obtain the intuitively self-consistent result 



F+ - 



1-IAo 
A^ 



Tf ) i?r Eq. 



(59) 



A similar mode-matching condition for the back reference 
plane and incident field can be found using the same pro- 
cedure. 



E. Michelson Interferometer 

The propagation schematic of a generalized Michelson 
interferometer (MI) is shown in Fig. ^. The back (AR) 
parallel input and output ports are coupled to either a 
mirror or FPI Mi through a propagation distance Lig, 
and the front (HR) perpendicular ports are coupled to 
either a mirror or FPI Mii through a propagation dis- 
tance 1/26- We can simplify later calculations of the op- 
tical performance of a power-recycled gravitational wave 
interferometer if we first determine a transfer matrix for 
the MI, and then establish a procedure for satisfying the 
dark port condition in the presence of perturbations. As 
in Section |ll D| , we develop the MI transfer matrix in a 
basis-independent fashion, and postpone the discussion 
of the choice of an unperturbed b asis s et until we calcu- 
late steady-state fields in Section IIIC . 



1. Transfer Matrix 

By comparing Fig. |^ and Fig. ||, we can connect the MI 
input and output fields of frequency component q with 
the front parallel and back perpendicular fields of the 
beamsplitter using the associations 



E^. 



F+W 



Fj 



As shown in Fig. H, the back parallel fields are coupled 
through reflection from optical component Mi, and the 
front perpendicular fields are coupled through reflection 
from Mil- Therefore, for frequency component q, we 
have 



F+^ 
9 



KeieiAiOg) E+K and (60a) 
Ke26i^i^<i)E;^. (60b) 



where 



KeieiAu;) 



G6iRiiAuj)Gie, and(61a) 



i2 Au tq2 



G62i?n(Ac^)G; 



26- 



(61b) 



In Eq. ( plaD , Gie is the Gouy operator describing the 
forward propagation step from the back parallel output 
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2. Idealized Simulation of the Dark Port Condition 

Interferometric gravitational-wave Michelson interfer- 
ometers typically operate with the beamsplitter posi- 
tioned so that the back perpendicular (antisymmetric) 
output port generates a dark fringe. In the model MI 
developed in the previous section, we can adjust the 
static position zq of the beamsplitter to ensure that 
an unperturbed, mode-matched fundamental field (i.e., 
{to, n, q} — {0, 0, 0}) incident on the beamsplitter's front 
parallel input port is cancelled at the back perpendicu- 
lar output port. Given our conventions for the phases of 
the amplitude reflection and transmission coefficients of 
a high-reflectance coating, we choose zq = — 7r/2\/2fcp^ 
and obtain from Eq. ( ^3|) 



^6^" = 

Rt' = 



R^^ = iree+'^^^'<' A^, (64a) 
-iTQtl e-' -^^^ Sl si Aq , and (64b) 
-ir^tl e-' Si sl . (64c) 



FIG. 9: Schematic diagram of the reflection and transmission 
operators of a Michelson interferometer, deflned by Eq. (^3|) . 
Note that the internal reflection operators i?f (Alj,) and 
i?j^(Ati;5) may each re prese nt either the reflection operator of 
a mirror, given by Eq. ([38a|) , or that of a Fabry- Perot interfer- 
ometer, given by Eq. ( p6a| ) . The labels z'^ and z^ denote the 
input and output reference planes, respectively, near the HR 
surface of the beamsplitter , while and C,^ are the input 
and output reference planes near the BS AR surface. 



port of the beamsplitter TWe to the input port of tUi, 
while Ggi is the corresponding operator describing prop- 
agation from the output port of Mi to the back parallel 
input port of Mq. If the length of the vacuum sepa- 
rating A^e from AA\ is Lig, then the single- pass propa- 
gation time is rig = Lw/c = rgi. In Eq. (|61b|) , there 
is an equivalent set of propagation operators and times, 
with 1 — > 2 and I — > II. The internal reflection operators 
i?]~(Att') and i?^(Aa;) may each represent eithe r the re- 
flection operator of a mirror, given by Eq. ( |38a), or that 
of a Fabry-Perot interferometer, given by Eq. (56a). 

By substituting Eq. ( |60| ) into Eq. (|4l|), we can imme- 
diately construct a transfer matrix similar to Eq. ( p6f) for 
the MI as a monolithic optical element. Using Fig. ^ as 
a guide, we find: 



T^iiAug) R 



MI (^^9) 
(Aw,)_ 



T+ 



(62) 



where the individual reflection and transmission opera- 
tors are defined by the matrix product 



T^i(Ac^) R 



'MI 
,^Ml(A^) 



MlV 

T+ 



(Ao;)] 

Ml(^w) 

KeieiAuj) 
Ke26{Aiu) 





^6 


Re 




iRe 


^6 J 



R,^ 



Rl 



(63) 



If a compensation plate is placed in the vacuum between 
the beamsplitter and A^n, then (depending upon its ori- 
entation angle) it can be represented as either a mirror or 
a beamsplitter substrate with two antireflecting coatings 
and no curvature, and the corresponding operator can 
be inserted in to th e perpendicular propagator _ftrg26(Aw) 
given by Eq. ( 61t ) . 

Suppose that the MI has been initialized using a partic- 
ular set of configuration parameters, and in the absence 
of a thermal lens in any substrate the dynamic microdis- 
placements of the beamsplitter and the mirrors are set 
to zero. If the static beamsplitter position has been set 
to satisfy the dark port condition for q — 0, then an 
unperturbed TEMqo field incident on the front parallel 
port produces a vanishingly small output at the back 
perpendicular port. Insofar as it is possible, we wish to 
develop an algorithm that maintains this dark port con- 
dition when the system has been perturbed. 

Consider an arbitrary field incident on the front 
input port of the MI. At the back output port, the trans- 
mitted field can be found by simplifying Eq. (^) for 
(7 = 0, giving 



E+ = [D^iAze) + D2{Aze)] F+ 



(65) 



where we choose to extract from Eq. (|6j) the explicit 
phase offset arising from a nonzero value of Azg, and 
write 

DiiAze) = e+'^''^''Rt^Keie{0)Ti\ anc(66a) 
D2{Aze) = e-'^''^'''T+^Ke2e{0)Re^. (66b) 

The corresponding power measured by a square-law de- 
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tector at the dark port is therefore (within a constant) 



1 



^1^0^ D\{Aze)D2{Aze)Fn+c.c. 
d\{0) D2{0)Fo 
X cos (2V2 k Azq — (jjj 

where (j) is the phase angle of —F^^^ ^i(O) -D2(0) F^ . In 
order to minimize P(^{Azq), we require that the condi- 
tion d Pq (Azq) / d Azq be satisfied, or 



(67) 



PRFPMI, we develop a simulation of an idealized servo- 
controlled length locking system that incorporat es the 
dark port condition algorithm described in Section [I E S . 
We then demonstrate the determination of a set of unper- 
turbed basis functions that can be used to describe the 
electromagnetic intracavity and extracavity fields, and 
we develop a numerical algorithm for the computation of 
the steady-state fields under the influence of both geo- 
metric and thermal perturbations. After a detailed com- 
parison of the predictions of this model with those of a 
fast Fourier transform code set in two important special 
cases, we evaluate the optical performance of the baseline 
LIGO design. 



Az6 = 



(68) 



We will incorporate this procedure into a general algo- 
rithm for simulating a servo-controlled resonator length- 
locking system for an entire gravitational-wave interfer- 
ometer in Section IIIB. 



III. 



POWER-RECYCLED FABRY-PEROT 
MICHELSON INTERFEROMETER 



In previous sections, we have constructed the trans- 
fer matrices of the components needed to build a func- 
tionally complete description of a power-recycled Fabry- 
Perot Michelson interferometer (PRFPMI]^ The PRF- 
PMI propagation schematic shown in Fig. H is a simpli- 
fied form of the FPI schematic displayed in Fig. H, and 
will allow us to determine the reduced transfer matrix 
with a minimum of calculation. In this case, the mir- 
ror A^i is replaced by the power-recycling mirror A^s, 
and the mirror A^3 is replaced by the Michelson interfer- 
ometer A^Mi- After determining the corresponding en- 
hancement, reflection, and transmission operators of the 



Transfer Matrix 



Consider the schematic representation of the enhance- 
ment, reflection, and transmission operators of the inter- 
ferometer shown in Fig. |lO|. If we express these operators 
using the transfer matrices of A^i and A^mi, respectively 

then the in- 



described in Section 11 B S and Section II E 1 



tracavity field enhancements can be represented as 



F,- = H+iAu;,)F+, and 
F+ = H+iAu;,)F+. 



(69a) 
(69b) 



where is the amplitude of the external sideband field 
q incident on A^s, F^^^ is the amplitude of frequency com- 
ponent q incident on reference plane z^, and Fq is the 



corresponding field incident on reference plane Zg. By 
comparison with Eq. (Q), the enhancement operators 



are 



H+{Au;) = [1 



„i 2 Auj T(35 



(70a) 
(70b) 



where Ges is the Gouy operator describing the forward 



Given the PRFPMI schematic shown in Fig. 



propagation step from to Zg, and G5g is the cor- the enhancement operators H^{A(jj) and Hq{Au 



responding operator describing propagation from Zg to 
z^. If the length of the vacuum separating A^s from 
A^Mi is ^65 1 then the single-pass propagation time is 
res = Lq5 I c = T56 . The front transmission and reflection 
operators of t he p ower-recy cling mirror, and i?^, are 
given by Eq. ( |37a| ) and Eq. (38a), respectively. The front 
reflection operator of the MI, i?j^j(Aijj), is given by the 
appropriate element of Eq. (p3) . 



obtain the output fields E^ and E'^ as 



E- 



RYpoi^uj,)F+ 



and 



(71a) 
(71b) 

where the transmission and reflection operators are 



K = 



T+o{Auj,)F: 



T+o(Ac^) = T+i(Au;)iJ6+(Aw), and 



' IFO 
^IFO 



i?T,.0(^^) 



Rt 



(Au;) 



(72a) 
(72b) 
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FIG. 10: Schematic diagram of the enhancement, reflection, and transmission operators of a power-recycled Fabry-Perot 
Michelson interferometer, defined by Eq. (TO) and Eq. , respectively. 



B. Idealized Simulation of Servo-Controlled 
Resonator Length Locking 

Although the transfer matrix schematics of the FPI 
and the PRFPMI — shown in Fig. | and Fig. |l^, re- 
spectively — are similar, the resonant structure of the 
MI and the dark port condition significantly complicate 
simulations of an idealized resonator length-locking sys- 
tem. Clearly, we need to adjust the position of the power 
recycling mirror AI5 to maintain a "round-trip" reso- 
nance condition, but the round trip must include both 
resonances in the Michelson arms and reflection from a 
beamsplitter that has been infinitesimally displaced to 
force compliance with the Michelson dark port condition. 
In fact, we again choose to compute a microdisplaccmcnt 
A^s for directly from a numerical determination of 
the net round-trip phase accumulated by the lowest-loss 
carrier eigenmode of the entire resonator. Beginning and 
ending at the reference plane zf, the perturbed round- 
trip propagator at the carrier frequency for the PRFPMI 
shown in Fig. is 



find both the eigenvalue Aq with the largest magni- 
tude and the corresponding lowest-loss eigenvector 
Eq . Set the Michelson input field = Ges Eq 
in Eq. (|67|), a nd the n use Eq. (|6|) and the discus- 



sion in Section [I E 2 to compute the new beamsplit- 
ter position Azg that satisfies the dark port condi- 
tion for Eq. Since Eq is generally not an eigen- 
vector of iiriFo(0, Azg), iterate this step to deter- 
mine an eigenvector Eq and a BS displacement Azg 
such that successive values of the largest-magnitude 
eigenvalue |Aq| differ by no more than a suitably 
small convergence threshold (typically 10^^^ in our 
simulations). 



3. As in the case of the FPI discussed in Section 11 D 2 



the eigenvector £'q is also an eigenvector of the ma- 
trix operator /'iriFo(Az5, Azg), with the eigenvalue 
As ^ \\'q\ e»(<^o-2feA25)^ Hence, if we choose 



Az5 



2fc' 



(74) 



i4:iFo(A2;5, Aze) = e 



—i 2 k Azs 



Tg "if6i6(0)T^ 



(73) 



! 2^2 A2 



'i?6 "^626(0) i?6 



Ggs i?5 



wher e w e have chosen to extract from Eq. (38a) and 



Eq. (64a) the explicit phase offsets arising from nonzero 
values of A25 and Azg, respectively. 

A close examination of Eq. ( [73| ) reveals that we can 
separate the idealized servo-locking algorithm into three 
consecutive stages, each with its own optimization pro- 
cess. 



1. Follow the procedure described in Section [I D 2 to 
adjust individually the positions of Ali and M.2 
to maximize the round-trip enhancement for the 
lowest-loss carrier eigenmode of each FPI. 

2. Select an initial value for Azg, and then compute 
the carrier eigenvalue spectrum of Kypo{Q, A zg) to 



then the net round-trip phase accumulated by Eq 
will be zero, and that mode will satisfy the reso- 
nance condition. 

As in the case of the FPI, in our simulations we maintain 
a record of the most recently chosen lowest-loss eigenvec- 
tor, and in the case of degeneracy we choose the eigen- 
mode having the largest value of the inner product with 
the previous eigenvector. 

Suppose that we have followed this procedure and 
computed the eigenvector Eq of the matrix operator 
-K^iFo(0c)/2 ^' ^-^e)' representing a maximally resonant 
field amplitude at reference plane in Fig. O. We can 
mode-match this field by choosing an input field that 
optimally couples to E q through the enhancement oper- 
ator given by Eq. (70a). If we solve the matrix equation 
Eq = H+{0) F+ for F+ we obtain 



F+ - 



1-IAol 
AS 



i?5 Eq, 



(75) 



in agreement with Eq. (p£ 



17 



C. Computation of Steady-State Fields 

In previous sections, we have been largely unconcerned 
with the basis chosen to represent the operators which 
comprise the transfer matrices of the optical elements of 
the gravitational-wave interferometer. Instead, we have 
focused on the details of specifying the operator algebra 
needed to describe the physics of mirror perturbations. 
In this section, we intend to define a computational algo- 
rithm which will allow us to determine the self-consistent 
steady-state fields everywhere in a realistic model of a 
thermally loaded interferometer. We begin by describing 
a procedure to define an unperturbed transverse spatial 
basis that can be used to represent all intracavity and ex- 
tracavity electromagnetic fields, and then we specify the 
initialization process we have used to prepare the com- 
putational model for perturbations. Finally, we present 
the simple iterative steady-state solution algorithm we 
have used to obtain convergence from the complete set of 
nonlinear coupled propagation equations describing the 
PRFPMI intracavity fields. 

In the Michelson interferometer, the beamsplitter cou- 
ples two optical systems (each consisting of a vacuum 
region followed by either a mirror or a Fabry-Perot in- 
terferometer) that may have significantly different mode- 
matching requirements relative to some common refer- 
ence plane. In principle, the determination of a set of 
eigenmodes that is common to both optical systems is 
a tedious numerical exercise, particularly in the context 
of the power-recycling scheme shown in Fig. |l|. How- 
ever, in practice, our simulations do not require the 
identification of an ab initio set of common eigenmodes. 
Rather, we seek a collection of self-consistent unper- 
turbed basis functions capable of accurately represent- 
ing the transverse spatial features of the laser field any- 
where in the gravitational-wave interferometer. As we 
show here, these basis functions can be eigenmodes of 
the unperturbed parallel FPI only, propagated unidirec- 
tionally throughout the entire system. 

Referring to Fig. |l|, Fig. ^, and Fig. ^, we begin by ig- 
noring all apertures everywhere in the system, and choos- 
ing the length and mirror curvatures of the unperturbed 
(e.g., infinite-aperture, spherical- mirror, and cold) par- 
allel FPI. These geometric configuration parameters de- 
fine a set of standing-wave eigenmodes with a particular 
spot size and a wavefront radius of curvature at i that 
matches that of Ali prior to thermal loading. We se- 
lect these eigenmodes as the spatial basis functions that 
will describe the transverse laser field everywhere in the 
system. 

Referring to Fig. |l| and Fig. |l^, we position another 
mirror at the location of the power recycling mirror 
(PRM) Ai^to couple optically the two Fabry-Pcrot inter- 
ferometers, and we propagate the lowest-loss (i.e., "fun- 
damental" ) basis function out of the parallel FPI through 
the beamsplitter to the intracavity input reference plane 
of the PRM. We calculate the value of the curvature 
1 / Rp of the extracted fundamental FPI eigenmode at , 



and set the unperturbed curvature of A^s to that value. 
We also compute the spot radius W5 of the fundamental 
basis function at z^ so that we can determine the per- 
turbed curvature mismatch matrices using the true 



curvature l/R^ in Eq. (jUj) and Eq. (|lj). 

Following the same unidirectional approach, we prop- 
agate the fundamental basis function from the PRM in- 
tracavity output reference plane z^ to the intracavity 
reference planes Z2 of A^2 and then z^ of A^4 in the 
perpendicular FPI. At each of these two reference planes, 
we compute the curvature and spot radius of the corre- 
sponding laser field, and then set these values as the ap- 
propriate unperturbed parameters of these mirrors. As 
in the case oi M5, any discrepancies between the un- 
perturbed curvatures of the propagated basis functions 
and the true curvatures of the perpendicular mirrors are 
then captured by the perturbation matrices defined 
by Eq. (0). 

We complete our specification of the unperturbed ge- 
ometry of the interferometer by propagating a selected 
set of unperturbed transverse spatial eigenmodes from 
the parallel FPI throughout the PRFPMI. This proce- 
dure defines a (possibly incomplete) set of spatial basis 
functions with which intracavity fields in the PRFPMI 
can be specified. Furthermore, by propagating the in- 
tracavity basis functions out through both the dark port 
and A^5, we can use the same unperturbed basis for the 
extracavity fields F"*", , and E~ shown in Fig. ^ In 
general, the input field will not be described solely by 
the outward-propagated fundamental eigenmode of the 
parallel FPI; instead, it will be represented by a linear 
combination of some subset of the modes available in the 
unperturbed basis. 

The spherical-mirror, infinite-aperture, cold approxi- 
mation of the unperturbed gravitational-wave interfer- 
ometer lends itself to a convenient representation of intra- 
cavity fields by power-orthogonal Hermite-Gauss modes. 
However, even when the mirrors remain spherical and 
thermal lensing is ignored, this basis cannot represent all 
possible perturbed fields with arbitrarily high accuracy 
under all circumstances. For example, in cases where 
either the recycling cavity is geometrically unstable or 
a thermal lens is so strong that the sign of the curva- 
ture of a wavefront propagating through the substrate 
changes, then we expect that the number of stable un- 
perturbed basis functions needed to represent that field 
will become prohibitively large. Furthermore, if either 
finite apertures or non-spherical mirrors are used to de- 
fine the unperturbed basis numerically, it is possible that 
residual diffraction of a non-Hermite-Gauss fundamen- 
tal basis function from Ai^ (followed by propagation to 
Aie) will generate a field profile that has a mean spot ra- 
dius that differs from the value computed using the field 
profile obtained by propagating the lowest-loss parallel 
eigenmode from A^i to Alg- In this biorthogonal 
basis can be chosen, but special care must be taken to en- 
sure that the eigenmode expansions converge, When a 
self-consistent basis can be chosen, and the correspond- 
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ing perturbation expansion of the physical observables 
under simulation converges, the matrix representation of 
the problem can allow computations which are orders of 
magnitude faster than those performed by FFT codes. 

Now that we have defined a self-consistent set of spa- 
tial basis functions at every reference plane in the LIGO 
PRFPMI, we can introduce the perturbations by setting 
the apertures and curvatures of the mirrors and beam- 
splitter to their true values, and then computing the ma- 
trices that characterize the effects of aperture diffraction, 
wavefront-mirror curvature mismatch, and thermal fo- 
cusing on laser fields expressed as linear combinations 
of th e un perturbed basis functions. Initially, we solve 
Eq. ( |69a|) in the limit to obtain the perturbed 

fields of the cold interferometer, and then we choose a 
relatively low but finite input power (e.g., 0.1 W) for the 
input field. The steady-state intracavity fields under this 
small thermal load can be found using a straightforward, 
self-consistent procedure: 

1. Using the current intracavity fields, compute the 
powers absorbed in the mirror and beamsplitter 
substrates and coatings, and update the substrate 
thermal OPD matrix elements given by Eq. ( |3^ ) 
and Eq. (|4^), which depend nonlinearly on these 
powers. 



Following the procedure outlined in Section pi B| , 
adjust the micropositions of A4i, M2, Me, and 
A^5 to achieve optimum resonant phase conditions 
for the lowest-order carrier mode in the each of the 
three optical cavities. 



3. Recompute the intracavity fields using Eq. ( 69a ) 



4. Repeat the previous steps until the recycled power 
i \ E^q\ has stabilized. 

Note that we do not propagate the field from one ref- 
erence plane to another throughout the interferometer. 
In this case, even when using the simulated phase con- 
trol system, the fields tend to converge after a number 
of iterations implied by the photon storage lifetimes of 
the arm cavities. Furthermore, we do not use more com- 
plicated nonlinear gradient-search methods to determine 
the intracavity fields because the number of variables 
(field basis function coefficients) is extremely large, and 
the variables at any single location in the interferome- 
ter is linked to those at all other locations by nontrivial 
round-trip propagators. Instead, the solution procedure 
we have detailed above can provide convergence to a few 
parts in 10^ in only a few iterations. 

The steady-state fields arising from higher input pow- 
ers can be found by initializing the convergence process 
with scaled field amplitudes computed at a lower power. 
For example, if a stable solution has been found at 2 W 
total input power, then the same collection of intracav- 
ity fields, scaled by a factor of two, would be reasonable 
initial values for an input power of 8 W. Similarly, if the 



TABLE IV: High-reflection (hr) and anti-reflection (ar) loss 
parameters for the initial LIGO optical elements used in our 
simulations. Here is the optical power absorption that con- 
tributes to the substrate thermal loss, and Sc is the scattering 
loss. 



Parameter 


Loss (ppm) 


aiir 


1.0 




1.0 


^hr 


90.0 




600.0 



TABLE V: Physical parameters for the initial LIGO optical 
elements used in our simulations. Each substrate is fused 
silica, with parameters given by Table | and a radius a = 
25 cm. Here the total loss l^r = <^hr + Shr, where a^r and s^r 
are the power absorption and scattering losses in the high- 
reflection coating, respectively. 



Mirror 


Rm (m) 


h (m) 


r^ 


t^ 


Ml, M2 


14 571.0 


0.10 


0.970000 


I - - Ihr 


M3, Mi 


7 400.0 


0.10 


1 — — Ihr 


1.5 X 10"^ 


Ms 


9 999.8 


0.10 


0.985020 


1 - - Ihr 


Ma 


00 


0.04 


0.499975 


1 - - Ihr 



input power is held constant, but some other parameter 
of the PRFPMI configuration (such as the PRM curva- 
ture) is varied, the intracavity fields can be re-used from 
one parameter value to the next. In this way, a set of so- 
lutions for a variety of different interferometer parameter 
sets can be constructed in a few minutes using a typical 
desktop computer. 

An object-oriented numerical computer model of the 
initial LIGO interferometer based on these algorithms 
has been implemented using the class mechanism of the 
MATLAB programming language. [2^ The resulting col- 
lection of MATLAB source code files, hereafter referred 
to as "Melody," is available publicly for examination and 
execution. [ETI 



D. Comparison with a Fast Fourier Transform 
Model 



In certain limited cases, the predictions made by 
Melody can be compared with those of a Fortran imple- 
mentation of a Fast Fourier Transform (FFT) model of 
the initial LIGO interferometer . ||l7| The FFT model was 
designed to investigate geometric requirements and tol- 
erances for initial LIGO optical components, and allows 
either a measured or simulated static OPD phase map 
to be included in each refiective surface. By con trast, 
the Melody solution method described in Section IIIC 



allows the evolving fields everywhere in the interferome- 
ter to alter the local phase maps using nonlinear thermal 
lensing, so that both the field and mirror are distorted 
by their reciprocal interaction. For the FFT the paraxial 
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approximation is used in order to calculate the propaga- 
tors between mirrors, as matrix operators in the spatial 
frequency domain. The lengths are optimized in order to 
achieve a stationary locked configuration, and the power 
stored inside the arms and the recycling cavity is eval- 
uated. The final fields obtained may be analyzed for 
modal composition. In contrast, the Melody program 
starts with a fixed number of paraxial modes. The ma- 
trix elements coupling the modes are then analyti cally 
calculated using the algorithms described in Section [I B 
and Section II C , with the thermal effects perturbatively 



evaluated at every iteration until stationarity is achieved. 

While the basic optics and wave-equation assumptions 
are the same, the two programs have markedly different 
implementations. In the case of the FFT code, the carrier 
and the sidebands are simulated by two separate, consec- 
utive runs so that the cavity lengths are optimized for the 
carrier, and the dimensions of the recycling cavity are op- 
timized for the sidebands. As discussed previously, the 
idealized Melody control system sets the cavity lengths 
by ensuring a round trip phase of zero for the lowest-loss 
eigenmode. The FFT solutions scale linearly with ini- 
tial field power, but in Melody, the power is shared by 
carrier and sidebands; the carrier and sideband fields are 
separately calculated, and their combined power at each 
iteration affects the thermal lens distortion nonlinearly. 
The FFT features fine-grained meshing of fixed phase 
map information, being intended primarily to study ar- 
bitrary mirror aberrations. On the other hand. Melody 
is intended to much more rapidly study coarser-grained 
effects with specific emphasis on the non-linear effects of 
thermal focusing. Significantly, the computational hard- 
ware requirements for the codes is dramatically different: 
a typical FFT run requires 30-60 minutes on a super- 
computer, while a typical Melody computation requires 
only a few seconds to stabilize. 

The parameters that we used in our simulations are 
shown in Table |^ and Table 0. The length of each arm 
cavity was chosen to be 3999.01 m, and the RF modu- 
lation frequency cji/27r = 24.493 MHz (with a modula- 
tion depth of 0.44) was adjusted slightly in the Melody 
computations to guarantee an antiresonance in the arm 
cavity to within machine precision. The current design 
value of the initial LIGO common length (the average of 
the distances between the power recycling mirror and the 
input test masses) is 9.188 m, and the differential length 
(or "Schnupp") asymmetry (the difference between the 
distances between the parallel and perpendicular ITM 
and the PRM) is Ls = 0.15 m. The distance between 
the beamsplitter and the PRM has the value 6.253 m, 
and all simulations were performed using the Nd:YAG 
laser wavelength Aq = 1.0642 fim. As a first comparison, 
both codes were used to simulate the unperturbed ini- 
tial LIGO interferometer, with perfect mode matching, 
infinite-diameter optics, and no thermal distortion. The 
agreement between computations of reflected and trans- 
mitted power with analytical calculations are within 0.2% 
for both the FFT code and Melody. 
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FIG. 11: A comparison between the predictions of computer 
code based on our model (Melody) and those of the FFT 
model in the case of simple aperturing. The power enhance- 
ment (or "gain") at the parallel/inline FPI reference plane 
z^" in Fig. H is plotted as a function of the aperture of Mg. 
The "clip approximation" represents the result expected by 
assuming that the loss arises from simple single-pass aperture 
clipping at Ms- 
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FIG. 12: A comparison between the predictions of Melody 
and those of the FFT model in the case of unstable curvature 
reduction. The gain at the recycling mirror reference plane 
At in Fig. is plotted as a function of the curvature of 
Ai2- Note that — in the absence of thermal loading — the 
recycling cavity becomes unstable for the resonant sidebands 
at a radius of curvature of 14480 m. 



In the absence of thermal loading, we have stud- 
ied the dependence of the parallel (inline) arm cavity 
power buildup on the diameter of the ETM Ms- A 
comparison between the predictions of Melody and the 
FFT code tests the suitability of the former's basis- 
function-expansion approach to calculating interferom- 
eter power levels, particularly when the cavity losses are 
large enough to require the inclusion of a large number of 
modes. In Fig. Ill], we plot the ratio of the total circulat- 
ing optical power at reference plane in Fig. |^ to the 
total carrier input power iji^g^p. When the 28 lowest- 
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loss Hermite-Gauss modes (corresponding to m + n < 6) 
of the parallel FPI are used, the agreement is excellent 
even at an ETM diameter of only 15 cm, which imposes 
a round-trip diffractive loss of approximately 2000 ppm. 
For comparison, we have included a trace displaying the 
stored power in the parallel arm expected when the intra- 
cavity loss is due to simple single-pass aperture clipping. 
Note that this naive assumption underestimates the true 
diffractive loss by a factor of about 2.14 for this beam 
spot size at A^3. 

In the thermally loaded simulations that are described 
in Section [II E , the increase in thermal focusing power 



of the ITM substrates improves the geometric stability 
of the power recycling cavity by increasing the appar- 
ent ITM radii of curvatures. However, if we reduce the 
radius of curvature (ROC) of either ITM below the de- 
sign value of 14.571 km, the recycling cavity becomes 
unstable, resulting in a large reduction in stored power. 
Previous studies of this phenomenon using the FFT code 
have shown that — for small changes in ITM curvature — 
only the sidebands suffer this power degradation; the car- 
rier resonance condition in the corresponding arm cavity 
serves to stabilize the recycling cavity mode by attenu- 
ating the anti-resonant higher-order modes. Therefore, 
since our model uses stable resonator eigenmodes as the 
spatial basis functions for the field expansions, this is 
a much more strenuous test than the aperture reduc- 
tion example, requiring a substantially greater number 
of transverse modes. In Fig. ^ we show the dependence 
of the power recycling cavity gain on the perpendicular 
(outline) ITM radius of curvature. The Melody results 
for the carrier agree essentially perfectly with those of the 
FFT code even when we employ only the 28 lowest-loss 
Hermite-Gauss modes of the parallel FPI. However, to 
obtain reasonably small discrepancies between the pre- 
dictions of the sideband behavior in the highly unsta- 
ble regime, we must include all basis functions satisfying 
m + n < 20, or 231 modes. 



E. Initial LIGO Optical Performance Simulations 

We now turn to simulations of the basic initial LIGO 
interferometer in the presence of thermal loading of the 
optics. Because of the high power levels maintained in 
the interferometer to suppress shot noise, the absorption 
of power in the ITM substrates and coatings and the sub- 
sequent thermal focusing must be incorporated into the 
interferometer design. Using Melody, we explore in detail 
the effects of this focusing on the optical performance of 
the interferometer. Since the recycling cavity is operat- 
ing in a region of geometric stability under full thermal 
load, we have run our simulations using the 28 lowest-loss 
Hermite-Gauss modes (corresponding to m -I- n < 6) of 
the parallel FPI. 

If the radius of curvature of the power recycling mir- 
ror A^5 is chosen to be the cold-cavity mode- matched 
value shown in Table the thermal distortions caused 



by the substrates in the recycling cavity will destroy this 
mode-matching and significantly reduce the available re- 
cycling gain for the sidebands. In Fig. |l^ and Fig. |l^, 
we plot the sideband recycled powers |F5,±ip given by 
Eq. ( |69a| ) as a function of the PRM ROC ^5 for the ini- 
tial LIGO TEMoo operating input laser power of 6.5 W, 
with 5.9 W delivered by the carrier when the sideband 
modulation depth is 0.44. The presence of the arm cavi- 
ties significantly reduces the sensitivity of the carrier to 
these perturbations, resulting in an approximately con- 
stant recycling gain of 30 when the parameters specified 
in Table |, Table 0, and Table |l^ are used. Note that 
the optimum mode-matching of the parallel and perpen- 
dicular segments of the recycling cavity occurs when the 
PRM ROC has increased to approximately 15 km. 

One of the most striking details of Fig. ^ and Fig. |lj 
is the significant difference in the recycling gain of the up- 
per {+A0J1) and lower (— Ao^i) sideband fields. Such a 
sideband imbalance can lead to significant gravitational- 
wave detection noise signatures that are qualitatively ab- 
sent when the balance is exact. Since the ratio Auji/luq 
is immeasurably small for RF sidebands and (by design) 
there is a very high degree of geometrical symmetry be- 
tween the Michelson branches of the interferometer, we 
naively expect the sideband fields to interact identically 
with geometrical distortions as long as these distortions 
are frequency-independent at the RF scale. In principle, 
this symmetry is broken only by the macroscopic value 
of the Schnupp asymmetry Ls = isi — ^52, where L^j is 
the optical path length separating the M5 and Mj. For 
gravitational- wave interferometers, Ls is much smaller 
than the Rayleigh length, so both sidebands should inter- 
act nearly identically with either arm, except for a phase 
accumulation proportional to the distortion-independent 
factor AujiLs- Using both Melody and the FFT code, we 
have determined that our simulations strictly show excel- 
lent sideband balance if either — > or the net branch 
distortions are identical. However, with Melody we have 
determined that the branch distortions cannot be identi- 
cal when the thermal lens in the beamsplitter substrate 
is significant; this source of additional distortion in the 
parallel arm requires a significant displacement of the 
beamsplitter to maintain the dark port condition, and 
is therefore the primary cause of the imbalance shown 
in Fig. ^ and Fig. |lj. In fact, as seen in Fig. |l^, the 
FFT code reproduces this behavior when the field-mirror 
curvature mismatch at the perpendicular ITM becomes 
large, causing a significant difference between the net dis- 
tortions of the two arms. Further studies using Melody 
have shown that it is possible to improve the sideband 
balance by supplementing the procedure outlined in Sec- 
tion IIIB with an additional step that readjusts the po- 



sitions of both A^s and A4q, at the expense of increasing 
the TEMqo carrier output power. 

For the purposes of our discussion here, we have cho- 
sen a value of the PRM ROC that minimizes both the 
fundamental carrier power and the sideband imbalance 
at the output port. This strategy produces the optimum 
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FIG. 13: Initial LIGO recycled powers as a function of the 
PRM radius of curvature. The total TEMoo input laser power 
is 6.5 W, with 5.9 W delivered by the carrier. 



FIG. 15: The optimum initial LIGO PRM radius of curvature 
as a function of total TEMoo input laser power. 
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FIG. 14: Initial LIGO dark port powers as a function of the 
PRM radius of curvature. The total TEMoo input laser power 
is 6.5 W, with 5.9 W delivered by the carrier. 



initial LIGO PRM radius of curvature as a function of to- 
tal TEMoo input laser power shown in Fig. ^ Again, we 
have assumed a modulation depth of 0.44, which pushes 
about 10% of the total available laser power into the side- 
bands. The corresponding optimized initial LIGO recy- 
cled and output sideband powers are shown in Fig. |l^ 
and Fig. |l^, respectively. Note that at each value of the 
input power in these two plots, the corresponding PRM 
ROC shown in Fig. |l5|has been used. Although the resid- 
ual imbalance of the recycled sideband fields remains at 
the 10% level at an input power of 10 W, the imbalance 
at the antisymmetric port is negligible. 

As an example of the operating characteristics of a 
fixed interferometer design point. Fig. |l^ and Fig. |l^ 
give the recycled and output powers at the carrier and 
sideband modulation frequencies for an interferometer 
optimized for a total laser input power of 6.5 W. The 
optimum PRM ROC is 15.075 km, and, as shown in 



Fig. 18, the recycling gain for the carrier is about 30. Fig- 



ure |19| predicts that, at 6.5 W, the sidebands are approxi- 
mately balanced, with about 130 mW in each field, while 
the total carrier power emerging from the dark port is 
about 180 /LtW. Approximately 90% of this carrier power 
is stored in higher-order modes. 

The behavior of the idealized PRFPMI phase-control 
system is illustrated in Fig. |o[ The PRM (TWs) is 
initially pulled in toward the beamsplitter to compen- 
sate for the longitudinal phase introduced by the mis- 
match between the curvatures of the field emerging from 
the parallel FPI and the PRM itself. However, as the 
strengths of the thermal lenses in the ITM substrates 
increase, the mirror is pushed away from the beamsplit- 
ter to maintain the carrier resonance. Similarly, as the 
astigmatic thermal lens develops in the beamsplitter sub- 
strate, the beamsplitter moves slightly toward the up- 
per left in Fig. |^. Indeed, as we see in Fig. even in 
the y plane the effective focal length of the beamsplit- 
ter substrate is six times larger than that of either ITM 
substrate; nevertheless, the resulting asymmetry between 
the two Michelson arms is numerically detected and com- 
pensated. 



IV. CONCLUSION 

We have developed highly detailed steady-state ana- 
lytical and numerical models of the physical phenomena 
which limit the intracavity power enhancement of recy- 
cled Michelson Fabry-Perot laser gravitational-wave de- 
tectors. We have included in our operator-based formal- 
ism the effects of nonlinear thermal lensing due to power 
absorption in both the substrates and coatings of the 
mirrors and beamsplitter, the effects of any mismatch 
between the curvatures of the laser wavefront and the 
mirror surface, and the diffraction by an aperture at each 
instance of reflection and transmission. We have taken 
great care to preserve numerically the nearly ideal longi- 
tudinal phase resonance conditions that would otherwise 
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FIG. 16: Optimized initial LIGO recycled sideband powers as 
a function of total TEMoo input laser power. At each value 
of the input power, the PRM radius of curvature shown in 
Fig. has been used. 



FIG. 19: Initial LIGO dark port powers as a function of total 
TEMoo input laser power. 




FIG. 17: Optimized initial LIGO output sideband powers as 
a function of total TEMoo input laser power. At each value 
of the input power, the PRM radius of curvature shown in 
Fig. has been used. 



FIG. 20: Initial LIGO PRM and BS microdisplacements as a 
function of total TEMoo input laser power. 
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FIG. 18: Initial LIGO recycled powers as a function of total 
TEMoo input laser power. The PRM radius of curvature is 
15075 m, the optimized value for an input power of 6.5 W. 



FIG. 21: Initial LIGO effective (i.e., quadratic) thermal focal 
lengths as a function of total TEMoo input laser power. 
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be provided by an external servo-locking control system. 
Once the mirrors and beamsplitters have been analyzed 
and their transfer matrices constructed, more compli- 
cated optical components and structures (such as Fabry- 
Perot and Michelson interferometers, and the LIGO de- 
tector) can be modeled using simple matrix composition. 

Our formalism relies on the validity of an expansion 
of the intracavity electromagnetic field at any reference 
plane in a necessarily incomplete set of possibly biorthog- 
onal basis functions that are derived from eigenmodes of 
one of the Fabry- Perot arm cavities. In two specific cases, 
we have checked the results produced by Melody — a 
fast MATLAB implementation of our model against 
those generated by an FFT program used to model opti- 
cal imperfections in LIGO. We have found that Melody 
can accurately predict the behavior of the carrier even 
in highly perturbed cold resonators using a modest num- 
ber of spatial basis functions in orders of magnitude less 
time than is required by the FFT code. The same is 
true for the sideband fields stored in the power recycling 
cavity, except in the case where the recycling cavity is sig- 
nificantly geometrically unstable. Even in this situation 
(which does not arise when the interferometer is ther- 
mally loaded), accuracy can be improved substantially 
by including more modes in the simulation. 

We have discovered and described in broad terms the 
conditions under which the power stored in the two side- 
band fields can become unbalanced, potentially affecting 
the noise sensitivity of the gravitational-wave interferom- 
eter in the detection band. At this point, we believe that 
an additional post-process step can be introduced into 



our simulated phase-locking algorithm that will allow a 
reduction of the sideband imbalance at the expense of a 
previously optimized optical property of the carrier field. 

In the future, we will include the effects of thcrmoelas- 
tic surface deformation in our model, and apply it to pro- 
posed configurations of the advanced LIGO interferom- 
eter. Other systematic perturbations and imperfections 
(e.g., mirror tilt and substrate inhomogeneities) can be 
included easily by incorporating the appropriate opera- 
tors into the transfer matrices describing reflection and 
transmission for the mirrors and beamsplitter. In addi- 
tion, we intend to develop algorithms for the simulation 
of the optical response of an interferometer to gravita- 
tional radiation, with the intention of estimating the de- 
tection sensitivity of advanced LIGO under full thermal 
load. We anticipate that the performance of the Melody 
program will not be significantly affected by these mod- 
ifications. 
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